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We have developed the capability to determine accurate harmonic spectra for multielectron atoms 
within time-dependent J?-matrix (TDRM) theory. Harmonic spectra can be calculated using the 
expectation value of the dipole length, velocity or acceleration operator. We assess the calculation 
of the harmonic spectrum from He irradiated by 390 nm laser light with intensities up to 4 x 10 14 
Wcm -2 using each form, including the influence of the multielectron basis used in the TDRM code. 
The spectra are consistent between the different forms, although the dipole acceleration calculation 
breaks down at lower harmonics. The results obtained from TDRM theory are compared with 
results from the HELIUM code finding good quantitative agreement between the methods. We find 
that bases which include pseudostates give the best comparison with the HELIUM code, but models 
comprising only physical orbitals also produce accurate results. 
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I. INTRODUCTION 

In recent years, harmonic generation (HG) has become 
one of the richest veins of research for atomic, molecular 
and optical physics. Not only has HG enabled ultrashort 
light pulse generation but it has also given rise to 
a series of very sensitive measurements of molecular 
atomic and even electronic dynamics [|[. The sensitive 
nature of HG has made it increasingly important to de- 
velop accurate methods of modelling the process. 

Many studies aimed at describing HG make use of the 
single active electron (SAE) model [1, @, a significant 
simplification which allows for the efficient computation 
of harmonic spectra. SAE methods have been used to 
probe the relationship between atomic structure and HG. 
For instance, the Cooper minimum in Argon has been 
linked with the minimum in the photoionization spec- 
trum caused by a zero dipole moment between the p 
ground state wavefunction and the d wavefunction of the 
photoionized electron for a photon energy of around 48 
eV @,0j. The minimum is observed to exist in both the 
photoionization spectrum and the harmonic spectrum, 
and is easily described by the SAE method as it does 
not depend on the interactions between different elec- 
trons. However, there have been various studies carried 
out in molecular systems where multielectron dynamics 
are found to be of great importance (HQ. Even in atomic 
systems there are features of photoionization, and hence 
harmonic spectra, that are the result of electronic interac- 
tions which require a multielectron description 0, 0-[H| ■ 

Over the last few years, we have developed time- 
dependent i?-matrix (TDRM) theory to model the in- 
teraction of atoms with short, intense laser pulses, main- 
taining a full description of the multielectron dynamics 
involved TDRM has recently been extended to 

account for harmonic generation, and this capability was 
demonstrated in showing how autoionizing resonances 
can affect the harmonic spectrum of Argon. The ap- 
pearance of the autoionizing resonances in these spectra 



is a consequence of multielectron dynamics: the inter- 
ference between the response of 3p and 3s electrons to 
the laser field @ . These calculations represent an impor- 
tant shift in thinking on HG: the multielectron nature of 
the process is reflected in the theoretical approach, and 
while there are many processes that can be adequately 
described using SAE methods, there are many for which 
this, more rigorous, description may be required. 

The determination of the harmonic spectrum can pro- 
ceed through the calculation of the time-dependent ex- 
pectation value of the dipole, dipole velocity or dipole ac- 
celeration operator. At present there is discussion about 
which of these operators offers the best prediction of the 
harmonic response for a single atom. Recent work [Infill ] 
has suggested that there is a natural connection with the 
dipole velocity, while, commonly, the dipole acceleration 
operator is used [r7Hl9j . especially for the description 
of high order harmonics as better resolution can be ob- 
tained for the high energy peaks [20f. Much early work 
in the field used the dipole length [2lT [22j , and up un- 
til this point the description of HG in TDRM has been 
restricted to using this operator @. We note that the 
use of these various forms has been verified only within 
the SAE approximation, and hence we seek herein to 
verify the independence of HG with respect to the use 
of the dipole, its velocity or acceleration in a multielec- 
tron system. We also assess which form offers the most 
numerically stable method, particularly when used with 
a limited multielectron basis set. Studies assessing the 
propagation of the wavefunction have demonstrated that 
to obtain the most accurate results for a limited basis in 
the TDRM approach the laser field is best described in 
the dipole length gauge [23]. On the other hand, time- 
propagation in SAE calculations is commonly performed 
by describing the laser field in the velocity gauge. This 
difference indicates that we cannot necessarily rely upon 
knowledge gained from SAE calculations for the assess- 
ment of TDRM calculations. 

We have extended the TDRM method to calculate the 
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harmonic spectrum from the dipole velocity and acceler- 
ation operators simultaneously with the dipole operator 
spectrum. In this paper we cover the major theoretical 
aspects of this extension, and apply the TDRM codes to 
He in a 390nm laser field. Helium is chosen for three 
reasons. Firstly, the simple structure allows for the sys- 
tematic varying of the multielectron basis functions, the 
impact of which has been assessed for TDRM in terms of 
photoionization [23[, but not for HG. Secondly, the ab- 
sence of a closed core simplifies the calculation of dipole 
acceleration matrix elements, and hence we can compare 
spectra in all three forms. Finally, using He allows us 
to benchmark our approach against a proven alternative 
method: we compare our results with those obtained us- 
ing the HELIUM code [H. 



II. THEORY 

A. TDRM Theory 

The TDRM approach is an ah initio nonperturbative 
theory for describing ultrafast atomic processes. Details 
of the method can be found in [12, HH , so we only give 
a short overview here. The time-dependent Schrddinger 
equation for an atom containing (N + 1) electrons is 

i^(X N+1 ,t) = H(t)V(X N+ll t). (1) 

The Hamiltonian, H, contains both the non-relativistic 
Hamiltonian of the N + 1-electron atom or ion in the ab- 
sence of the laser field and the laser interaction term. The 
laser field is described using the dipole approximation in 
the length form, and is assumed to be linearly polarized 
and spatially homogeneous. This form provides the most 
reliable ionization yields when only a limited amount of 
atomic structure is included [23j. 

We propagate a solution of the time-dependent 
Schrddinger equation ^ on a discrete time scale with time 
step At in a Crank-Nicolson scheme. We can write the 
wavefunction at a time t q+ i in terms of the wavefunction 
at the previous time step t q : 



-{H m +E)* t 



(2) 



Here the imaginary energy E is defined as 2i/At and H m 
is the Hamiltonian at the midpoint of the time interval, 

tq+l/2- 

In i?-matrix theory, configuration space is partitioned 
into an inner and outer region. In the inner region, all 
electrons are within a distance a- ln of the nucleus, and 
full account is taken of all interactions between all elec- 
trons. In the outer region, an ionized electron moves be- 
yond the boundary a- ln , and thus exchange interactions 
between this electron and the electrons remaining close 
to the nucleus can be neglected. The ionized electron 
then moves in only the long-range multipole potential of 
the residual A-electron core and the laser field. 



Following 13[ we can evaluate Eq. ([2j at the boundary 
£ii n as a matrix equation 



F (a in ) = R(a in )F (a in ) + T (a in ) . 



(3) 



in which the wavefunction F, at the boundary is de- 
scribed in terms of its derivative, F, plus an inhomo- 
geneous vector, T, arising from the right hand side of 
Eq. @. The R- matrix, R, connects the inner and outer 
region wavefunction at the boundary ai n . 

Given an inner region wavefunction, R and T are eval- 
uated at the boundary a- m . Subsequently, they are prop- 
agated outwards in space up to a boundary, a out where 
it can be assumed that the wavefunction F has vanished. 
The wavefunction vector F is set to zero and then propa- 
gated inwards to the inner region boundary. Once F has 
been determined at each boundary point, the full wave- 
function can be extracted from the i?-matrix equations. 
We can then iterate the procedure using Eq. ((2J. 



B. Harmonic generation 

The electric field produced by an accelerating charge 
is given, using the non-relativistic Lienard-Wiechert po- 
tentials in the far field limit, by 



E(t) = k(*l>(t) 



[p z ,H] 



ih 



i/){t))+ keE User (t), (4) 



where e is the electronic charge, z is the laser polarization 
axis, A: is a proportionality constant, p z is the canonical 
momentum and i?i as0 r is the electric field of the laser 
pulse. We can write 



[Pz-H] 



ih 



dt 



m)\pz\m), (5) 



and it follows that 



E(t) oc d(t) 



dt 2 



(m\*w))- 



(6) 



The power spectrum of the emitted radiation is then 
given, up to a proportionality constant, by |d(cj)| 2 - the 
Fourier transform of d(t) squared. 

Although the radiation produced is proportional to the 
dipole acceleration, it is common practice in HG calcula- 
tions to calculate d(u>), i.e. to use the expectation value 
of the dipole length instead. This is because a simple re- 
lationship exists between d and d which can be extended 
to include the dipole velocity form: 



|dH| 



|d(o;)| 2 = |d(w) 



(J) 



Therefore the harmonic response of a single atom can 
be expressed in terms of the expectation value of the 
dipole operator 



d (*) = (*(*) I -«*!*(*)>> 



(8) 
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or of its velocity 



d(t) = -(*(t)|-ez|*(t)>, 



or acceleration 



d(*) = ^(*(t) 



!*(*)>> 



(9) 



(10) 



where z is the total position operator along the laser po- 
larization axis. 

As discussed in [23| the TDRM code can use either 
the length or velocity gauge for the propagation of the 
wavefunction. While, in keeping with the findings of (23[, 
we use the length gauge, we can still utilize the dipole 
velocity matrix elements produced by the i?-matrix suite 
of codes which 'seed' the TDRM code. Thus we can store 
both z and dz/dt and use Eqns. © and © directly for 
the determination of the time- varying expectation values 
of the dipole operator and the dipole velocity. 

However, in order to calculate the expectation value of 
the dipole acceleration we cannot use Eq. (fTO)) directly. 
Instead, using Ehrenfest's theorem, it is possible to write 
the dipole acceleration as 

•; , x ,dH. ,eZcos9, „ T . . ._, , , 
d (t) = (-^r) = <-^-> - eN elec (V\E(t)\V), (11) 

where Z is the nuclear charge, r the total position opera- 
tor, 9 the angle between f and z and N e i ec the number of 
electrons. The second term in Eq. (JTTJ is often seen with- 
out this factor of N e i ec as in the SAE approximation it is 
just 1. We can make a small change to the way the radial 
integrals are calculated in the i?-matrix suite which al- 
lows the calculation of (1/r-r) instead of (r) . Then we can 
use Eq. (fTTj) to calculate the dipole acceleration. Thus, 
we can now simultaneously calculate harmonic spectra 
using the dipole length, velocity and acceleration opera- 
tors. The propagation of the wavefunction is still carried 
out in the length gauge. 

The use of the acceleration form will however be re- 
stricted to He like targets. The use of Ehrenfest's theo- 
rem (in Eqs. ([5]) and (jlip) requires that the wavefunction 
be exact, or close to it. For general multielectron systems 
we normally impose a fixed core where (at least) the first 
two electrons are restricted to a single orbital. Imposing 
this restriction means that the electronic repulsion is not 
fully described. More precisely, if the orbital of electron 
ei is fixed, and the orbital of electron e2 is not, then the 
action on e2 will not necessarily equal minus the reaction 
on e±. Thus, the commutator 



(Pi +P2), ■ 



1 



|ri - r 2 

may not be guaranteed to vanish. On the other hand 

1 



(12) 



will still vanish. Thus, while the expectation 
value ([rcos#, H]) can be calculated accurately, 
([[r cos ,H}) can not, rendering Ehrenfest's theorem 
untenable. Thus, the comparisons we employ for the 
simple He test case which follows can be extended to 
general multielectron systems only for the dipole length 
and velocity forms. 



(ri+r 2 ), 



ri - r 2 



(13) 



C. Calculation parameters 



The one-electron basis used for describing the residual 
He + in the inner region consists of orbitals expressed in 
terms of B splines. The residual He + ion is represented 
through a series of models of increasing complexity [23| . 
The basic model consists of only the He + Is state, which 
we call IT (1 True state). We also use two models com- 
prising six states. The first is built using true orbitals Is, 
2s, 2p, 3s, 3p, 3d (6T) and the other using 5 pseudo or- 
bitals and the true Is orbital: Is, 2s, 2p, 3s, 3p, 3c? called 
6P, (6 states with Pseudostates) . Pseudostate models 
have been found to be more accurate in the time propa- 
gation of the He wavefunction responding to short light 
fields, especially in the velocity-gauge description of the 
light field. Pseudostate models may thus provide a better 
basis for the description of the ionization and HG pro- 
cesses, provided that these processes are not affected by 
artificial resonances introduced by the pseudostates. 

The inner region radius is set at 20 a.u. which is suf- 
ficiently large to contain the residual ion for each model 
we use. The outer region boundary is set at 600 a.u. to 
prevent any reflections of the wavefunction for the dura- 
tion of the short laser pulse employed. The set of contin- 
uum orbitals contains 80 B splines for each angular mo- 
mentum, £, of the continuum electron up to a maximum 
value, L max = 19. Convergence testing was carried out 
retaining angular momenta up to a value of L max = 27 
and, while changes in the harmonic spectra are observed, 
they occur at energies beyond the cutoff- outside the re- 
gion of interest here. The outer region is divided into 
sectors of 2 a.u. containing 35 9th order B-splines per 
channel. The time step used in the wavefunction propa- 
gation is 0.1 a.u. 

We use 390 nm laser pulses, consisting of a 3 cycle 
sin 2 ramp-on followed by two cycles at peak intensity, 
followed by a 3 cycle sin 2 ramp-off (3-2-3). We also cal- 
culate spectra for different pulse shapes and find that 
while the spectra change, the comparisons between them 
are generally described by the results presented below for 
the 3-2-3 pulse. There is one important exception to this 
general observation, which is discussed in Sec. MI Dl 
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FIG. 1: (Color online) The harmonic spectrum (up to a con- 
stant of proportionality) as calculated from the dipole accel- 
eration for He in a 390 nm, 4 x 10 14 Wcm -2 , 3-2-3 laser 
field, using as a model residual ion description, the Is state 
(black, dotted line), the Is, 2s, 2p, 3s, 3p, 3d states (red, 
dashed line), and the Is, 2s, 2p, 3s, 3p, 3d pseudostates (blue, 
solid line). The single state model provides a reasonable ap- 
proximation to the more detailed descriptions beyond the first 
harmonic where there is a large discrepancy between the spec- 
tra. 



FIG. 2: (Color online) The length form harmonic spectrum 
of He in a 390 nm, 4 x 10 14 Wcm^ 2 , 3-2-3 laser field, using 
as a residual ion description, the IT (black, dotted line), 6T 
(red, dashed line), and 6P (blue, solid line) models (See Fig. [TJ 
for details). Results from the IT model provide a reasonable 
approximation to results from the more detailed descriptions 
especially in the cutoff region. There is also good agreement 
for the first harmonic when compared with the large discrep- 
ancy in Fig. [T] 



III. RESULTS 
A. Comparison of various target states 

The harmonic response, as calculated from the expec- 
tation value of the dipole acceleration, of a He target 
in the IT, 6T and 6P configurations is shown in Fig. 
[TJ The spectra display the expected form- a pronounced 
first harmonic peak followed by a plateau of peaks at 
odd multiples of the fundamental photon energy, which 
decay exponentially beyond a cutoff. The cutoff of the 
plateau appears at a photon energy of approximately 45 
eV. The standard formula for the cutoff energy, I p +3.2U P 
[25j , where I p is the ionization potential and U p the pon- 
deromotive energy, is not necessarily appropriate in this 
wavelength and intensity regime. Nevertheless, for the 
current parameters, it predicts a cutoff energy of 42 eV. 
The observed cutoff is therefore not inconsistent with the 
cutoff formula. 

We can compare the spectra to assess how the descrip- 
tion of atomic structure affects the calculated HG spec- 
tra. The IT calculations are in better agreement with the 
more detailed calculations at higher energies, especially 
in the cutoff region between the 13th and 19th harmon- 
ics where agreement is within 30%. In the low energy 
region, and especially in the first harmonic, the spectra 
differ significantly- the first harmonic response in the IT 
model is 60 times greater than that in the 6P model. 
The inconsistencies in the lower harmonics between the 
IT and the 6T and 6P models imply that the low en- 
ergy harmonics in the dipole acceleration calculation are 
highly sensitive to changes in the atomic structure, and 



that in the higher energy cutoff region the details of the 
atomic structure are not as important. There is a factor 3 
difference in the first harmonic peak between the 6T and 
6P models. As pseudostates may better represent the 
changes to the ground state than true states, this differ- 
ence implies that the first harmonic is especially sensitive 
to the description of the ground state. The 6P spectrum 
shows a double peak structure at the 9th harmonic stage 
which the two true state models do not. 

As the first term in the dipole acceleration is pro- 
portional to 1/r 2 it is most sensitive to changes in the 
wavefunction at small r. If the description of the atomic 
structure close to the nucleus is not exact, this can lead 
to significant inaccuracies in the low energy region of the 
spectra calculated from the dipole acceleration, especially 
the first harmonic peak. Figure [5] shows the same har- 
monic spectra as Fig. Q]but in this case the spectra are 
calculated from the dipole length operator. In this form 
the harmonics are far less sensitive to the details of the 
atomic structure close to the nucleus as can be seen by 
the excellent agreement between the three spectra at the 
first and third harmonics (within 20%). 

In fact the agreement between the spectra from the dif- 
ferent target states is generally better in the length and 
velocity forms than in the acceleration- except for the 
9th and 11th harmonics, the agreement between the 6T 
and 6P dipole- length spectra is within 20%. This further 
highlights that the dipole acceleration is especially sen- 
sitive to the description of atomic structure. The main 
difference between the three spectra appears again in the 
9th harmonic. This difference is very similar to the dif- 
ference seen in Fig. [JJ in which the dipole acceleration 
was used to determine the harmonic spectrum. This in- 
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FIG. 3: (Color online) The harmonic spectrum of a pseu- 
dostate (6P) He target in a 390 nm, 4 x 10 14 Won" 2 3-2-3 
laser field, as calculated from the dipole length (black, dotted 
line), velocity (red, dashed line) and acceleration (blue, solid 
line). Agreement to within 20% is found between all three 
spectra up to the 19th harmonic peak (60 eV). The spectra 
diverge beyond this. 

dicates that this difference originates from the different 
bases used, rather than the choice of operator for the de- 
termination of the harmonic spectrum. This topic will 
be discussed further in section MI CI 



B. Comparison of dipole length, dipole velocity 
and dipole acceleration forms 

As has been addressed in the previous section, TDRM 
theory can calculate harmonic spectra from the dipole 
length, dipole velocity or dipole acceleration operators. 
We have already seen how the dipole acceleration is sen- 
sitive to the description of the atomic structure, particu- 
larly when it comes to the low energy region of the spec- 
trum. 

Figure |3] shows the harmonic spectrum of 6P He in a 
390 nm, 4 x 10 14 Wcm" 2 laser field as calculated using 
the dipole length, velocity and acceleration forms of the 
dipole matrix elements. The pseudostates model gives a 
more accurate description of the changes in the ground 
state due to the laser pulse, and hence should give a 
more accurate picture of the harmonic spectrum than 
the true state model. In terms of the agreement between 
the spectra this holds true, as the 6P model gives a con- 
sistent agreement between the three different approaches 
to calculate the harmonic spectrum where the 6T model 
breaks down at low harmonics. For the 6P model the 
three spectra agree within 20% at every harmonic peak 
up to the 19th, well into the cutoff region. In the 6T 
there is agreement within 20% between the dipole length 
and velocity spectra, but the dipole acceleration spec- 
trum differs by 60%, 30% and 40% in the first, third and 
fifth harmonics respectively. 

Regardless of which model is used, the three spectra di- 
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FIG. 4: (Color online) The high energy harmonic spectrum of 
true state (6T) He target in a 390 nm, 4 x 10 14 Wcm -2 3-2-3 
laser field, as calculated from the dipole length (black, dotted 
line), velocity (red, dashed line) and acceleration (blue, solid 
line). The leftmost harmonic shown is the 19th, above which 
the spectra diverge. 

verge beyond the 19th harmonic (Fig. [4|, with the dipole 
length spectrum becoming noisy and the dipole acceler- 
ation spectrum displaying a few more weak harmonics 
decaying into noise. The dipole velocity spectrum on the 
other hand displays a second plateau of peaks not seen 
in the other spectra. These peaks are not predicted clas- 
sically, and their absence from the other spectra implies 
that they are spurious. This implies that the length and 
velocity forms are reliable, but only in an energy range 
up to and including the cutoff region. This is especially 
important as for general multielectron targets the accel- 
eration form will be prohibitively sensitive to the limita- 
tions in the description atomic structure (See Sec. Ill CI) . 
However, by using both the dipole length and dipole ve- 
locity operators it is possible to obtain reliable harmonic 
spectra for multielectron systems using the TDRM ap- 
proach. 



C. Comparison with HELIUM 

Having demonstrated that the TDRM method is self 
consistent within a certain energy range, we now seek to 
benchmark our results against those from a proven alter- 
native method. The HELIUM method [2J] uses direct 
numerical integration of the full-dimensional TDSE to 
describe a two-electron system. By solving the TDSE 
directly, no significant approximations are made, and 
thus all important multielectron effects are included. 
This makes HELIUM an excellent code against which 
to benchmark TDRM. 

Figure [5] shows the length form harmonic spectra pro- 
duced by the 6T and 6P models of He alongside that 
produced by the HELIUM code, for a target in a 390 
nm, 4 x 10 14 Wcm -2 , 3-2-3 laser field. At the harmonic 
peaks the agreement is very good. The 6P and HELIUM 
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FIG. 5: (Color online) The harmonic spectra as calculated 
from the dipole length produced from the 6T (black, dotted 
line) and 6P (red, dashed line) He models for TDRM, and 
from the HELIUM code (blue, solid line). The 6P and 6T 
spectra agree with the HELIUM spectrum to within 20% and 
30% respectively up to the 21st harmonic peak (except at the 
9th and 11th harmonics for the 6T spectrum). Inset: Both 
the 6P and HELIUM models have a structured peak at the 
9th harmonic. The IT (green circles) and 6T spectra do not. 



spectra agree to within 20% up to the 21st harmonics, 
while the 6T spectrum is within 30% except at the 9th 
and 11th harmonics 

The inset in Fig. [5] shows detail in the 9th harmonic 
from the three calculations, and the TDRM IT model. 
The 6P model and HELIUM spectra show a structured 
peak which the IT and 6T do not. The ponderomotivc 
energy in the laser field shifts the He ground state down 
by around 5.7 eV, shifting the ls3p bound state into the 
vicinity of the 9th harmonic peak. The presence of a 
bound state has been shown to give rise to such struc- 
ture in the harmonic peaks [|| . It is useful to notice that 
the 6T model may not describe the changes to the He + 
ground state in the laser field as accurately. Thus, the 
shift of the ls3p state peak may differ and consequently 
we do not observe the double peak structure in the 9th 
harmonic for the 6T spectrum. The IT model does not 
account for any changes to the He + ground state, and 
differences between the IT model and the other mod- 
els are thus even larger. Expansion of the basis set in 
the TDRM approach thus leads to a harmonic spectrum 
which gets closer to the benchmark harmonic spectrum 
obtained using the HELIUM code. 

The agreement for the TDRM velocity form spectrum 
is even better: within 15% when comparing the velocity 
form, 6P, TDRM spectrum with the length form HE- 
LIUM spectrum. The excellent agreement between the 
spectra serves to give weight to the results obtained from 
both methods. The sensitivity of the harmonic spectra 
to the description of atomic structure makes it even more 
remarkable that the two methods overlap, especially in 
the low energy region. Figure [6] shows the low energy 
region of the velocity form spectrum obtained from the 
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FIG. 6: (Color online) The low energy region of the dipole 
velocity harmonic spectrum produced from the 6P He model 
in TDRM (black, dotted line), and the length form harmonic 
spectra from the HELIUM code (red, dashed line), and its 
SAE derivative (blue, solid line). The TDRM and HELIUM 
spectra are indistinguishable, but the SAE spectrum overes- 
timates the harmonic spectrum at low energies. 

6P model TDRM code, alongside the length form spec- 
tra from the HELIUM code and from an SAE simplifi- 
cation derived from the HELIUM code (26[. The three 
spectra agree well in the first harmonic, whereas the ac- 
celeration form spectra (not shown) vary widely. This 
confirms that the dipole velocity and length are signif- 
icantly less sensitive to the description of atomic struc- 
ture close to the nucleus, and are probably more reli- 
able in the low energy, especially first harmonic, region. 
The gridspacing in HELIUM and the limited basis set in 
TDRM impose constraints on the calculations very close 
to the nucleus. These constraints make it likely that the 
acceleration form spectra are less reliable in the first har- 
monic, which could give rise to the discrepancy between 
the two methods. 

In the third and fifth harmonics the SAE model 
markedly overestimates the harmonic spectra obtained 
from both the TDRM and HELIUM models, which are 
in excellent agreement with each other. This implies that 
the SAE model is not sufficient to describe low energy 
harmonic spectra, and that the lowest energy harmon- 
ics are significantly more sensitive to atomic structure. 
We note that in the plateau and cutoff regions the SAE 
spectrum is in good agreement with the full HELIUM 
spectrum, lending justification to the use of the SAE 
approximation for investigating the generation of higher 
harmonics in He. 



D. Comparison of various pulse lengths 

To probe the effect of the laser pulse profile on the har- 
monic spectra, as well as the 3-2-3 (3 cycles, sin 2 ramp 
on, 2 cycles peak intensity, 3 cycles sin 2 ramp off) pro- 
file, we ran calculations for various longer pulses, namely 
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FIG. 7: (Color online) The dipole velocity harmonic spectra 
produced from the 6P He model for a 3-2-3 (blue, solid line) (3 
cycle, sin 2 ramp on, 2 cycles peak intensity, 3 cycle, sin 2 ramp 
off) and a 5-4-5 pulse (red, dashed line). Both pulses have a 
peak intensity of 4 x 10 14 Wcm -2 and a wavelength of 390 
nm. The longer 5-4-5 pulse gives rise to narrower harmonic 
peaks, but the peak values are still within 20% of the 3-2-3 
spectrum. Inset: There is a significant difference between the 
two spectra in the 9th harmonic. 



5-2-5, 3-4-3 and 5-4-5 pulses. Broadly speaking, while 
the spectra themselves change (with narrowing peaks for 
the longer pulses), the comparisons between the IT, 6T 
and 6P models, with the HELIUM results or between 
the dipole length, velocity and acceleration forms do not 
change significantly. Figure[7]shows the spectra produced 
by a 3-2-3 and a 5-4-5 laser pulse. The peak values are 
not changed significantly by the different pulse profile, 
but the longer 5-4-5 pulse gives rise to narrower peaks 
and greater contrast. This gives a greater energy reso- 
lution between different peaks. Therefore the broad 9th 
harmonic peak in the 3-2-3 spectrum in Fig. [7] is fur- 
ther broadened by the presence of the nearby ls3p bound 
state, whereas the narrower peak arising from the 5-4-5 
pulse is isolated from any nearby atomic structure. 



IV. CONCLUSIONS 

Wc have extended the calculation of harmonic spectra 
in TDRM theory by determining these spectra through 
the time-varying expectation value of the dipole length, 
dipole velocity and dipole acceleration operators, and ap- 
plied the adapted codes to He irradiated by a 390 nm, 
4 x 10 14 Wcm~ 2 laser field. We have compared the 
spectra calculated using each form, assessed the effect of 
changing the multielectron basis set used to describe the 
residual ion, and benchmarked our results against those 
obtained from the HELIUM method. 

We have shown that for harmonic photon energies up 
to and including the cutoff region the TDRM method 
provides results which are both self-consistent (between 
dipole length, velocity and acceleration forms) and con- 



sistent with an independent approach. The favorable 
comparison between the TDRM and HELIUM methods 
in the velocity and length form spectra implies that the 
present approach can provide excellent results. Care 
must be taken in the lower harmonics especially if us- 
ing the dipole acceleration operator where the sensitivity 
to inaccuracies in the description of the atomic structure 
can seriously affect the reliability of the spectra obtained. 

For general multielectron systems we can perform the 
calculations using both the dipole length and velocity, 
and compare the two spectra in order to establish bounds 
on the reliability of the results. Both methods give ex- 
cellent agreement for He well into the cutoff region. The 
divergence of the spectra beyond this occurs at energies 
which arc usually outside the region of interest. 

We have also probed the advantages of the various 
residual ion descriptions, which can be used within the 
TDRM method, finding that smaller basis sets, such as 
the IT single target state, provide an efficient way of test- 
ing the code, and a reasonable approximation to the har- 
monic spectrum, but larger basis sets give more detailed 
spectra, as would be expected from their better descrip- 
tion of the atomic structure involved. We also find that 
the inclusion of pseudostates in the He + basis seems to 
lead to more accurate harmonic spectra. This is particu- 
larly noticeable when compared with the highly accurate 
HELIUM method. This is largely due to the more precise 
way in which the pseudostate model describes the vari- 
ations in the He + ground state in response to the laser 
field. 

However, the use of pseudostates for general multielec- 
tron atoms can be problematic. By introducing non- 
physical thresholds into the system, pseudo-resonances 
can show up in the harmonic spectrum. These inadver- 
tant features do not appear in the He case presented here, 
as the energies at which they become important are out- 
side the harmonic region of interest. For general mul- 
tielectron atoms, this is not necessarily the case. This 
does not mean that accurate calculations are not possi- 
ble for larger atoms. Pseudostates can be used as long as 
care is taken- with knowledge of the position of pseudo- 
thresholds, unphysical resonances can be identified and 
disregarded. Secondly, although the 6T He model is not 
as close to the HELIUM spectrum as the 6P, it is still 
within 30% at every harmonic peak except the 9th and 
11th (in the dipole length spectra). Physical orbitals can 
thus also be used to improve accuracy of harmonic spec- 
tra. The number of physical orbitals required may be 
larger than if pseudo-orbitals are used, but this is not 
a fundamental problem: it affects only the scale of the 
calculations. With careful analysis of, and comparison 
between, pure physical orbital and pseudostate models 
we can reliably assess the accuracy of harmonic spectra 
for general multielectron systems using TDRM theory. 

Furthermore, even models which use only physical or- 
bitals already offer significant gains over SAE models. A 
simple example of this is HG in Ar + . Harmonics pro- 
duced by Ar + ions have been suggested to be the source 
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of the highest harmonics observed from a neutral Ar tar- 
get [13, 28]. The presence of three low-lying, 3s 2 3p 4 Ar 2 + 
thresholds can have a significant effect on the harmonic 
spectrum, and hence interactions between channels as- 
sociated with these thresholds must be accounted for. 
These interactions are neglected in an SAE calculation, 
but would be accounted for in a TDRM calculation in- 
volving purely physical orbitals. 

We find that the reliability of the results is not signifi- 
cantly affected by the particular laser pulse profile used. 
We compared results for four different laser pulse profiles, 
finding that while the harmonic spectra differed between 
cases, the changes were consistent between the various 
target state models, and with the HELIUM code results. 
In cases where atomic structure gives rise to structure in 
the harmonic spectrum the laser pulse length may affect 
the way in which this is observed. The greater energy res- 
olution afforded by longer pulses can isolate the separate 
effects of atomic structure. 

The results presented are also consistent with those 
from various peak intensities. We calculated harmonic 
spectra for intensities between I x I0 14 Wcm~ 2 and 
4 x 10 14 Wcm -2 , finding that the results are largely con- 
sistent. At lower intensities the plateau region is severely 
truncated and so it is difficult to compare between the 
various spectra, but the agreement is still evident in the 
cutoff region. 

The TDRM method has been rigorously tested up to 



intensities of 4 x I0 14 Wcm -2 and at wavelengths up to 
390 nm, but requires a significant amount of develop- 
ment to extend beyond these limits. It will be interesting 
to compare these findings with those that will be deter- 
mined using the new RMT (i?-matrix with time) codes 
[29l . [30l | which may be better suited to address higher 
intensities and longer wavelengths. While the TDRM 
method has been proven to provide interesting insight 
into the multielectron nature of HG, it has thus far only 
been implemented for general multielectron atoms us- 
ing the dipole length operator Q. The next stage will 
be to apply TDRM at high intensities to systems other 
than He. While the dipole acceleration is too sensitive 
to the description of atomic structure to accurately de- 
scribe such atoms, the length and velocity forms are sta- 
ble enough to provide good results for general targets. 
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